Emergence of quasi-units in the one dimensional Zhang model 
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We study the Zhang model of sandpile on a one dimensional chain of length L, where a random 
amount of energy is added at a randomly chosen site at each time step. We show that in spite of 
this randomness in the input energy, the probability distribution function of energy at a site in the 
steady state is sharply peaked, and the width of the peak decreases as L _1//2 for large L. We discuss 
how the energy added at one time is distributed among different sites by topplings with time. We 
relate this distribution to the time-dependent probability distribution of the position of a marked 
grain in the one dimensional Abelian model with discrete heights. We argue that in the large L 
limit, the variance of energy at site x has a scaling form L~ 1 g(x/L), where g(£) varies as log(l/f) 
for small £, which agrees very well with the results from numerical simulations. 

PACS numbers: 64.60.av, 05.65+b 
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I. INTRODUCTION 

After the pioneering work of Bak, Tang and Wiesen- 
feld (BTW) in 1987 1], many different models for self- 
organized criticality have been studied in different con- 
texts; for review see [2-5]. Of these, models in the general 
class known as Abelian distributed processors have been 
studied a lot, as they share an Abelian property that 
makes their theoretical study simpler The original 
sandpile model of Bak et al. the Eulerian walkers 
model @, and the Manna model Q are all members of 
this class. Models which do not have the Abelian prop- 
erty have been studied mostly by numerical simulations. 
In this paper, we discuss the Zhang model [7|, which does 
not have the Abelian property. 

In the Zhang model, the amount of energy added at a 
randomly chosen site at each time step is not fixed, but 
random. In spite of this, the model in one dimension 
has the remarkable property that the energy at a site in 
the steady state has a very sharply peaked distribution 
in which the width of the peak is much less than the 
spread in the input amount per time step, and the width 
decreases with increasing system size L. This behavior 
was noticed by Zhang using numerical simulations in one 
and two dimension [3], and he called it the 'emergence of 
quasi-units' in the steady state of the model. He argued 
that for large systems, the behavior would be same as 
in the discrete model. Recently, A. Fey et al. Q have 
proved that in one dimension, the variance of energy does 
go to zero as the length of the chain L goes to infinity, 
but they did not study how fast it decreases with L. 

In this paper, we study this emergence of 'quasi-units' 
in one dimensional Zhang sandpile by looking at how the 
added energy is redistributed among different sites in the 
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avalanche process. We show that the distribution func- 
tion of the fraction of added energy at a site x' reaching a 
site x after t time steps following the addition is exactly 
equal to the probability distribution that a marked grain 
in the one-dimensional height type BTW model added at 
site x' reaches site x in time t. The latter problem has 
been studied recently [l0( ■ We use this to show that the 
variance of energy asymptotically vanishes as 1/L. We 
also discuss the spatial dependence of the variance along 
the system length. In the large L limit, the variance 
at site x has a scaling form L~ 1 g(x/L). We determine 
an approximate form of the scaling function which 
agrees very well with the results of our numerical simu- 
lations. 

There have been other studies of the Zhang model ear- 
lier. Blanchard et al. [ll[ have studied the steady state 
of the model, and found that the distribution of energies 
even for the two site problem is very complicated, and 
has a multi-fractal character. In two dimensions, the dis- 
tribution of energy seems to sharpen for lar ger L, but 
the rate of decrease of the width is very slow [12] • Most 
other studies have dealt with the question as to whether 
the critical exponents of the avalanche distribution in 
this model are the same as in the discrete Abelian model 
[Til [3] • A. Fey et al.'s results imply that the asymptotic 
behavior of the avalanche distribution in one dimension is 
indeed the same as in the discrete case, but the situation 
in higher dimension remains unclear [TIL HE] ■ 

The plan of the paper is as follows. In Section II, we 
define the model precisely. In Section III, we show that 
the calculation of the way the energy added at a site is 
distributed among different sites by toppling is same as 
the calculation of the time-dependent probability distri- 
bution of the position of a marked grain in the discrete 
Abelian sandpile model. This correspondence is used in 
Section IV to determine the qualitative dependence of the 
variance of the energy variable at a site on its position 
x, and on the system size L. We propose a simple ex- 
trapolation form that incorporates this dependence. We 
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check our theoretical arguments with numerical simula- 
tions in Section V. Section VI contains a summary and 
concluding remarks. A detailed calculation of the solu- 
tion of an equation, required in Section IV, is added as 
an Appendix. 



II. DEFINITION AND PRELIMINARIES 

We consider our model on a linear chain of size L. The 
sites are labelled by integers 1 to L and a real continuous 
energy variable is assigned to each site. Let E(x, t) be 
the energy variable at site x at the end of the time-step t. 
We define a threshold energy value E c , same for each site, 
such that sites with E(x, t) > E c are called unstable, and 
those with E(x, t) < E c are called stable. Starting from 
a configuration where all sites are stable, the dynamics 
is defined as follows. 

(i) The system is driven by adding a random amount of 
energy at the beginning of every time-step at a randomly 
chosen site. Let the amount of energy added at time t be 
A t . We will assume that all A's are independent, identi- 
cally distributed random variables, each picked randomly 
from an uniform interval 1 — e < A t < 1 + e. Let the site 
of addition chosen at time t be denoted by at- 

(ii) We make a list of all sites whose energy exceeds or 
becomes equal to the critical value E c . All these sites are 
relaxed in parallel by topplings. In a toppling, the energy 
of the site is equally distributed to its two neighbors and 
the energy at that site is reset to zero. If there is toppling 
at a boundary site, half of the energy at that site before 
toppling is lost. 

(iii) We iterate Step (ii) until all topplings stop. This 
completes one time step. 

This is the slow driving limit, and we assume that all 
avalanche activity stops before the next addition event. 
In this limit, the model is characterized by two parame- 
ters e and E c . In the limit e = 0, and 1 < E c < 2, the 
model reduces to the discrete case, where the behavior 
is well understood [l7|]. For non-zero but small e, the 
behavior does not depend on the precise value of E c . In 
fact, starting with a recurrent configuration of the pile, 
and adding energy at some chosen site, we get exactly 
the same sequence of topplings for a range of values of 
£ c ||. To be precise, for any fixed initial configuration, 
and fixed driving sequence (of sites chosen for addition 
of energy), whether a site x topples at time t or not is 
independent of E c , so long as we have 1 + e < E c < 2 — 2e. 
In the following, we assume for simplicity that E c = 3/2, 
and < e < 1/4. 

It was shown in [9|] that in this case, the stationary 
state has at most one site with energy E(x, t) — and all 
other sites have energy in the range 1 — e < E{x 1 1) < 1+e. 
The position of the empty site is equally distributed 
among all the lattice points. There are also some re- 
current configurations in which all sites have energy 
E(x,t) > 1 — e. In such cases, we shall say that the site 
with zero energy is the site L + 1. Then, in the steady 



state, there is exactly one site with energy equal to 0, 
and the L + 1 different positions of the site are equally 
likely. 

If E c does not satisfy the inequality 1 + e < E c < 
2 — 2e, this simple characterization of the steady state 
is no longer valid. However, our treatment can be easily 
extended to those cases. Since the qualitative behavior 
of the model is the same in all cases, we restrict ourselves 
to the simplest case here. 

It is easy to see that the toppling rules are in general 
not Abelian. For example, start with a two site model in 
configuration (1.6,2.0) and E c = 1.5. The final configu- 
ration would be (1.4, 0), or (0, 1.3), depending on whether 
the first or the second site is toppled initially. In our 
model, using the parallel update rule, the final configu- 
ration would be (1.0,0.8). A. Fey et al. ^ have shown 
that only in one dimension, for 1 + e < E c , the Zhang 
model has a restricted Abelian character, namely, that 
the final state does not depend on the order of topplings 
within an avalanche. However, topplings in two different 
avalanches do not commute. 



III. THE PROPAGATOR, AND ITS RELATION 
TO THE DISCRETE ABELIAN MODEL 

It is useful to look at the Zhang model as a perturba- 
tion about the e = limit. For sufficiently small e, given 
the site of addition and initial configuration, the top- 
pling sequence is independent of e. It is also independent 
of the amount of energy of addition A t , and is same as 
the model with e = 0, which is the 1-dimensional Abelian 
sandpile model with integer heights (hereafter referred to 
simply as ASM, without further qualifiers). We decom- 
pose the energy variables as 

E{x, t) = Niat[E(x, t)\ + erj(x, t), (1) 

where Nint refers to the nearest integer value. Then the 
integer part of the energy evolves as in the ASM. We 
write 

A t = 1 + eu t , for all t. (2) 

Here Ut is uniformly distributed in the interval [—1, +1]. 
The linearity of energy transfer in toppling implies that 
the evolution of the variables rj{x, t) is independent of 
e. Thus, r/(x,t) is a linear function of u t ; the precise 
function depends on the sequence of topplings that took 
place. These are determined by the sequence of addition 
sites {at} up to the time i, and the initial configuration 
Co- These together will be called the evolution history 
of the system up to time t, and denoted by Ht- We 
assume that at the starting time t = 0, the variables 
t](x, t — 0) are zero for all x, and the initial configuration 
is a recurrent configuration Cq of the ASM. Then, from 
the linearity of the toppling rules, we can write 7y(x, t) as 
a linear function of for 1 < t' < t, and we can write 
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for a given history Tit, 

t 

rj(x,t\{u t },Ht) = ^ G(x,t\a t ,,t',Ht)u tl . (3) 
t'=i 

This defines the matrix elements G(x, t\at>,t', Tit)- These 
can be understood in terms of the probability distribu- 
tion of the position of a marked grain in the ASM as 
follows. Consider the motion of a marked grain in the 
one dimensional height type BTW model. We start with 
configuration Co and add grains at sites according to the 
sequence {at}. All grains are identical except the one 
added at time t' , which is marked. In each toppling, the 
marked grain jumps to one of its two neighbors with equal 
probability. Consider the probability that the marked 
grain will be found at site x after a sequence of relax- 
ation processes at time t. We denote this probability as 
Prob(ir, t\a t i , t', Tit)- From the toppling rules in both the 
models, it is easy to see that 

G{x,t\a v ,t',H t ) =Proh(x,t\a t ,,t',H t ). (4) 

Averaging over different histories Tit, we get the proba- 
bility that a marked grain added at x' = a t ' at time t' 
is found at a position x at time t > t' in the steady 
state of the ASM. Denoting the latter probability by 
ProbASu(x,t\x' ,t'), we get 



G(x,t\x> =a t ,,t>,H t ) =Prob ASM (x,t\x',t'), (5) 

where the over bar denotes averaging over different his- 
tories TLt, consistent with the specified constraints. Here, 
the constraint is that TL t must satisfy a t > — x' . At other 
places, the constraints may be different, and will be spec- 
ified if not clear from the context. 

We shall denote the variance of a random variable £ by 
Uar[£]. From the definition in Eq. (JTJ) , it is easy to show 
that 

Var[E{x, t)} = L/{L + l) 2 + e 2 Var[r)(x, t)]. (6) 

Different Ut are independent random variables, also in- 
dependent of TC t and have zero mean. Let T^ar[w t ] = a 2 . 
For the case when u t has a uniform distribution between 
— 1 and +1, we have a 2 = 1/3. Then, from Eq. ©, we 
get 

t 

Var[r){x, t)} = a 2 ^ G 2 (x, t\at> , t', Ut). (7) 
t'=i 

As t — > oo, the system tends to a steady state, and the 
average in the right hand side of Eq. ((7|) becomes a func- 
tion of t — t' . Also, for a given t', all values of ay are 
equally likely. We define 

F(x,t) = - lim V G 2 (x,t< + T\x>,t',H t ). (8) 

Li t' — »-00 * ' 

x' 



Then, for large L, in the steady state (t large), the vari- 
ance of energy at site x is 1/L + e 2 T l 2 (x) 1 where 

oo 

E 2 (x) = lim VarMx, t)] = a 2 V F(x, r). (9) 

t—>oo ^ — ' 

T = Q 

We define S 2 to be the average of S 2 (a;) over x. 

^ = iE s2 ^)- ( 10 ) 

X 

Evaluation of G(x, t\x', t' ,Ttt) for a given history Tit and 
averaging over Ht is quite tedious for t > 1 or 2. For G, 
the problem has been studied in the context of residence 
times of grains in sand piles, and some exact results are 
known in specific cases [10(. For G 2 , the calculations are 
much more difficult. However, some simplifications occur 
in large L limit. We discuss these in the next section. 

IV. CALCULATION OF E 2 (x) IN LARGE-L 
LIMIT 

In order to find the quantity F(x,t) in Eq. (8), we 
have to average G 2 (x,t\x' ,t' ,Ht) over all possible histo- 
ries Tit, which is quite difficult to evaluate exactly. How- 
ever, we can determine the leading behavior of F(x, r) in 
this limit. 

We use the fact that the_path of a marked grain in 
the ASM is a random walk [ljj| . Consider a particle that 
starts away from the boundaries, at x' = £L, with L 
large, and < £ < 1. If it undergoes r(TL t ) topplings 
between the time t' and t = t' + r under some particular 
history Ti. t , then its probability distribution is approxi- 
mately a Gaussian, centered at x' with width ^/r. Then, 
we have 

Using this approximation for G, summing over x' , we get 
Y^G 2 (x,t\x',t',n t )^ 1 (12) 

Thus, we have to calculate the average of 1/ y/r(H. t ) over 
different histories. Here r(TL t ) was defined as the num- 
ber of topplings undergone by the marked grain. Differ- 
ent possible trajectories of a marked grain, for a given 
history, do not have the same number of topplings. How- 
ever, if the typical displacement of the grain is much 
smaller than its distance from the end, differences be- 
tween these are small, and can be neglected. There are 
typically O(L) topplings per grain per avalanche in the 
model, and a grain moves a typical distance of 0[\[V) 
in one avalanche. Then, we can approximate r(Ti. t ) by 
N(x'), the number of topplings at x' . 

Let the number of topplings at x' at time steps r = 
0, 1, 2, . . . be denoted by N , Nx,N 2 , ■ ■ ■■ Then, N(x') = 
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No + N% + N% + • • • . It can be shown that the num- 
ber of topplings in different avalanches in the one dimen- 
sional ASM are nearly uncorrelated (In fact the correla- 
tion function between JVj and JVj varies as (1/L)\ l ~i\ .). 
By the central limit theorem for sum of weakly correlated 
random variables, the mean value of N grows linearly 
with r, but the standard deviation increases only as y/r. 
Then, for t»0, the distribution is sharply peaked about 
the mean, and ~ l/y/JW). 

Clearly, for r ^> 0, (N) — rn(x'), where n(x') is the 
mean number of topplings per avalanche at x' in the 
ASM, given by 



l(x = £L) = L£(l - 0/2. 



(13) 



The upper limit on r for the validity of the above ar- 
gument comes from the requirement that the width of 
the Gaussian be much less than the distance from the 
boundary, (without any loss of generality, we can assume 
that £ < 1/2, so that it is the left boundary ), else we 
cannot neglect eve nts wh ere the marked grain leaves the 
pile. This gives y/rn(x) <C £L, or equivalently, r <C £L. 
Thus we get, 

F{x,t)~ ^[ri£(l-0r 1/2 , for « r « £L, (14) 

where C\ is some constant. 

Also, we know that for r>L, the probability that the 
grain stays in the pile decays exponentially as exp(— t/L) 
[Tq| . Thus, G, and also G 2 will decay exponentially with 
r, for t>L. Thus, we have, for some constants C2 and 



F(x,t) 



C 2 



exp(— clt/L), for r » L. (15) 



It only remains to determine the behavior of Fix, r), for 
£L <C t <C L. In this case, in the ASM, there is a 
significant probability that the marked grain leaves the 
pile from the end. This results in a faster decay of G, and 
hence of F with time. We argue below that the behavior 
of the function F(x,t) is given by 



for £L « t « L, 



(16) 



where C3 is some constant. This can be seen as follows: 
Let us consider the special case when the particle starts at 
a site close to the boundary. Then n{x) is approximately 
a linear function of x for small x. Its spatial variation 
cannot be neglected, and Eq. (TH?)) is no longer valid. We 
will now argue that in this case 



G(x,f + T\x',t') 



1 -2 

X T 



exp(-x/r), 



for 0«t«I. 

(17) 

The time evolution of ProbASM(x, t\x', t') in Eq. ([5]) is 
well described as a diffusion with diffusion coefficient pro- 
portional to n(x) which is the mean number of topplings 
per avalanche at x in the ASM [ljj. For understanding 
the long-time survival probability in this problem, we can 



equivalently consider the problem in a continuous-time 
version: consider a random walk on a half line where 
sites are labelled by positive integers, and the jump rate 
out of a site x is proportional to x. A particle starts at 
site x — xo at time t — 0. If Pj(t) is the probability that 
the particle is at j at time t, then the equations for the 
time-evolution of Pj(t) are, for all j > 0, 

jP^t) = (i + l)P i+1 (i) + (i-l)P,_ 1 (t)-2jP J -(i). (18) 

The long time solution starting with Pj(0) = 5j >Xo is 

P,(t)~x t- 2 exp(-j/t) (19) 

for t 3> xo and large j. The probability that the particle 

survives till time t decreases as \ jt for large t. We have 

discussed the calculation in the Appendix. 

Using Eq. (O, we see that G(j, t' + t\xq, V) scales as 

Xq/t 2 . It seems reasonable to assume that G 2 will scale 
2 

as G . Then, each term in the summation for F[x, r) 
in Eq. ([5]) scales as Xq/t 4 , and there are r such terms, 
as the sum over Xq has an upper cutoff proportional to 
t, and so F(x,t) varies as 1/t for L ^> t ^> Xq. This 
concludes the argument. 

We can put these three limiting behaviors into a single 
functional form that interpolates between these, as 



F(x,r) 



1 K exp(— ar/L) 



L T + By/TLZ(1-Z)- 



(20) 



where K, a and B are some constants. In Section V, 
we will see that results from numerical simulation are 
consistent with this phenomenological expression. 

Using this interpolation form in Eq. @, and con- 
verting the sum over r to an integration over a variable 
u = t/L, we can write 

> r s f°° , K exx)(— au) ,„ . 

E 2 (a; = £Z)~— / du , ' ■ 21) 

L Jo u + By/u£(l-S) 

This integral can be simplified by a change of variable 
au = z 2 , giving 



H z {x = ^L) 



(22) 



where K, B' are constants, and I(y) is a function denned 
by 



I(y) = 2 dz 
Jo 



exp(— z 2 ) 
z + y 



(23) 



It is easy to verify that I(y) diverges as log(l/y) for small 
y. In particular, we note that the exponential term in the 
integral expression for I(y) has a significant contribution 
only for z near 1. We may approximate this by dropping 
the exponential factor, and changing the upper limit of 
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FIG. 1: (Color online) Scaling collapse of the probability dis- 
tribution Vl(E) of energy per site in the steady state for 
different systems of size 200, 500 and 1000. The distribution 
is well described by a Gaussian of width 0.136. 
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FIG. 2: (Color online) Scaling collapse of T, 2 (x)/a 2 at site x 
for systems of different length L. 



the integral to 1. The resulting integral is easily done, 



giving 



K'a 2 
L 



log 1 



BV£(i-0 



(24) 



where K' is some constant. Averaging E 2 (cc) over x, we 
get a behavior E 2 (x) ~ 1/L. Of course, the answer is 
not exact, and one could have constructed other inter- 
polation forms that have the same asymptotic behavior. 
We will see in the next Section that results from numer- 
ical simulations for T, 2 (x) can be fitted very well to the 
phenomenological expression in Eq. 



V. NUMERICAL RESULTS 
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FIG. 3: (Color online) The same plot in Fig. [2]resolved more 
at the left boundary of the model and taking x axis in log 
scale. 



We have tested our non-rigorous theoretical arguments 
against results obtained from numerical simulations. In 
Fig. [T] we have plotted the probability distribution 
Vl{E) of energy at a site, averaged over all sites. We 
used L = 200, 500 and 1000, and averaged over 10 8 dif- 
ferent configurations in the steady state. We plot the 
scaled distribution function Vl(E)/VL versus the scaled 
energy (E — E)*jL. A good collapse is seen, which veri- 
fies the fact that the width of the peak varies as L -1 / 2 . 

The dependence of the variance of E(x, t) on x is plot- 
ted in Fig. [2] for systems of length 200, 300 and 400. 
The data was obtained by averaging over 10 8 avalanches. 
We plot (L + A)£ 2 (x)/cr 2 versus x e ff/L e ff, where x e ff 
differs from x by an amount 6 to take into account the 
corrections due to end effects. Then, for consistency, L 
is replaced by L e ft = L + 26. For the specific choice of 
A = 5 ± 1 and 6 = 1.0 ± 0.2, we get a good collapse of the 
curves for different L. We also show a fit to the proposed 



interpolation form in Eq. (|2"3)). with K' = 1.00 ± 0.01 
and B' = 1.5 ± 0.2. We see that the fit is very good. 

In order to check the logarithmic dependence of £ 2 (x) 
on x for small x, we re-plot the data in Fig. [3] using 
logarithmic scale for x. We get a good collapse of the 
data for different L, supporting our proposed dependence 
in Eq. (El). 



VI. CONCLUDING REMARKS 

To summarize, we have studied the emergence of quasi- 
units in the one-dimensional Zhang sandpile model. The 
variance of energy variables in the steady state is gov- 
erned by the balance between two competing processes. 
The randomness in the drive i.e., the energy of addition, 
tends to increase the variance in time. On the other hand, 
the topplings of energy variables tend to equalize the ex- 
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cess energy by distributing it to the nearby sites. There 
are on an average 0(L 2 ) topplings per avalanche. Hence, 
in one dimension there are, on an average, O(L) topplings 
per site per avalanche. For large system size, the second 
process dominates over the first and the variance becomes 
low. We have shown that the variance vanishes as 1/L 
with increasing system size and the probability distribu- 
tion of energy concentrates around a non-random value 
which depends on the energy of addition. We have also 
proposed a functional form for the spatial dependence of 
variance of energy which incorporates the correct limit- 
ing behaviors, and matches very well with the numerical 
data. 

An interesting question is whether one can extend 
these arguments to the two-dimensional Zhang model. 
In this case, there are several peaks in the distribution of 
energies at a site, but there are some numerical evidences 
for the sharpening of the peaks as the system size is in- 
creased. However, as the number of topplings per site 
varies only as logL, the width is expected to decrease 
much more slowly with L, and the fluctuation effects can 
be much stronger. This remains an open question for 
further study. 



APPENDIX 

Here we discuss the solution of the Eq. (fT5| for the 
starting values given by 



P 3 it = 0) = o?-\ 



(A.l) 



We start with an ansatz Pj (i) = bt exp(— a t j), where both 
at and bt are functions only of t. This form satisfies the 
Eq. (fT8|) for all j, t > 0, if a t and b t satisfy 



da t 
dt 



2-e 



at _ „— a t 



(A.2) 
(A.3) 



To solve the Eq. |A~2|) . we first make a change of variable 
z = e~ at . In terms of z, the equation becomes dz/dt — 
(1 — z) 2 , which can be easily integrated to give 



t + A-1 
t + A ' 



(A.4) 



where A is an integration constant. To satisfy the initial 
condition in Eq. (|A.1[) . we choose 



A = (I -a)' 1 . 



Similarly, to solve the equation for bt, we use the form of 
e -a* given in Eq. (A4) and get 



db t l-2(t + A) 

—r- — bt- 



dt (t + A)(t + A - 1)' 
This can be integrated to give 

B 



(A.6) 



(t + A)(t + A-iy 



(A.7) 



where B is an integration constant. Then the probability 
can be written as 



{t + A-iy- 1 
(t + Ay+ 1 ' 



(A. 



To satisfy the initial condition at t = 0, we choose the 
integration constant B = (1 — a)~ 2 . Then, with these 
values of A and B, we have the solution for all j, t > 0, 
given by 



3-1 



W) = 



Now, as 4>j(a,t) satisfies the Eq. ([18 



[(1 - a)t + a] 
[(1 - a)t + 
4>j(a,t), say. 



(A.9) 



(n-1)! 9a"- 1 



(A.10) 



will also satisfy the equation for any natural number n. 
In addition, 



V>i,n(a = 0, t = 0) = Sj >n . 



(A.ll) 



Hence, we see that the solution of the Eq. (TTB)) . starting 
with Pj (t) = Sj^ n at t = is 



^•(*)=^j,n(a = 0,*) 



1 <9"-Vj(a,f) 



a=0; 



da 71 - 1 

(A.12) 

for all j, t > 0, where 4>j(a,t) is given in Eq. (|A.9|) and 
n is any natural number. 

It can be shown that for large t and j, the solution 



(A. 5) asymptotically becomes Pj(t) — nt exp(—j/t). 
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